%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Structural Model for Gallenstein 2024 JRI %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Results for Equality Treatment %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Table A.2 Model 1 %%
clc
clear
%D=xlsread('C:\Users\galle\Dropbox\Research\Working Papers\Inequality and Fairness Views\Data Analysis\Structural Model\EqualityOfferC.xls');
D=xlsread('C:\Users\galle\Dropbox\Research\Working Papers\Inequality and Fairness Views\Manuscript\Solo 2 Submissions\JRI Submission\Final Submission\Gallenstein_JRI_InequalityImpedesRiskManegement_Code&Data\EqualityOfferC.xls');
% initial guesses
mu    =    2.4611;    
sigma =    1.7077;    
Le    =    0.8720;    
a     =    0.9692;  
% for inequality treatments
v = [mu, sigma, Le, a]';
% Estimate parameters and Hessian Matrix
param = qnewton(@Loglikelihood_logit_EQ,v,[],D);
%Generating standard errors from the negative inverse Hessian matrix
hess = fdhess(@Loglikelihood_logit_EQ,param,D);
se = sqrt(-diag(inv(hess)));
[param se]
%% Calculating individual fairness view probabilities
% This section will produce fairness view predictions and output that data
% in excel format which then can be uploaded into STATA and used to
% generate table 5, 6, A.4, and A.5. After running this, use the STATA .do
% file: InequalityImpedeRiskManagement_Analysis.do 
mu = param(1);
sigma = param(2);
lambda_e= param(3);
a       = param(4);
w = D(:,1);
t1 = D(:,2);
t2 = D(:,3);
p = 0.7;
Pss = p*p;
Psf = p*(1-p);
Pfs = p*(1-p);
Y = 3.*w;
Yp= 3.*(12-w);
% Egalitarian - Equity
Fe1 = Y./2;
Fe2 = Yp./2;
% Egalitarian - Proportionate
Fp1 = (w/12).*Y;
Fp2 = (w/12).*Yp;
% Libertarian Proportionate 
Fl1 = (10/12).*Y;
Fl2 = (2/12).*Yp;

L_E=@(b) (1/(b*sigma*sqrt(2*pi)))*exp(-((log(b)-mu)^2)/(2*sigma^2)).*...
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*(t1.*Y).^(1-a)+Pfs*(1/(1-a)).*(t2.*Yp).^(1-a) - b.*(1./(2*Y)).*((t1.*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((t2.*Yp)-Fe2).^2)./...  
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)));

L_L=@(b) (1/(b*sigma*sqrt(2*pi)))*exp(-((log(b)-mu)^2)/(2*sigma^2)).*...
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*(t1.*Y).^(1-a)+Pfs*(1/(1-a)).*(t2.*Yp).^(1-a) - b.*(1./(2*Y)).*((t1.*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((t2.*Yp)-Fl2).^2)./...  
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)));


L = [integral(L_E,0,inf,'ArrayValued',true), integral(L_L,0,inf,'ArrayValued',true)];
Li  = lambda_e*L(:,1) + (1-lambda_e)*L(:,2);

likelihood = sum(log(Li))

FV_e_prob =  (lambda_e.*L(:,1))./Li;
FV_l_prob =  ((1-lambda_e).*L(:,2))./Li;
FV_prob=[L(:,1) L(:,2) FV_e_prob FV_l_prob D(:,4)]; 

T_EQ=table(FV_prob);
filename= 'EQ_probdata_CRRA.xlsx';
writetable(T_EQ,filename,'Sheet',1,'Range','D1')


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Results for Inequality Treatment %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Table A.2 Model 2 %%
clc
clear
D=xlsread('C:\Users\galle\Dropbox\Research\Working Papers\Inequality and Fairness Views\Data Analysis\Structural Model\InequalityOfferC.xls');
% initial guesses
mu    = 1.1135;    
sigma = 0.0051;    
Le    = 0.5966;    
Ll    = 0.2949;    
a     = 0.9681;   
v = [mu, sigma, Le, Ll, a]';
% Estimate parameters and Hessian Matrix
param = qnewton(@Loglikelihood_riskgame_logit_EvL2,v,[],D);
%Generating standard errors from the negative inverse Hessian matrix
hess = fdhess(@Loglikelihood_riskgame_logit_EvL2,param,D);
se = sqrt(-diag(inv(hess)));
[param se]
%% Calculating individual fairness view probabilities
% This section will produce fairness view predictions 
mu = param(1);
sigma = param(2);
lambda_e= param(3);
lambda_l= param(4);
a       = param(5);

w = D(:,1);
t1 = D(:,2);
t2 = D(:,3);

p = 0.7;
Pss = p*p;
Psf = p*(1-p);
Pfs = p*(1-p);
Y = 3.*w;
Yp= 3.*(12-w);

% Egalitarian - Equity
Fe1 = Y./2;
Fe2 = Yp./2;
% Egalitarian - Proportionate
Fp1 = (w/12).*Y;
Fp2 = (w/12).*Yp;
% Libertarian Proportionate 
Fl1 = (10/12).*Y;
Fl2 = (2/12).*Yp;

L_E=@(b) (1/(b*sigma*sqrt(2*pi)))*exp(-((log(b)-mu)^2)/(2*sigma^2)).*...
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*(t1.*Y).^(1-a)+Pfs*(1/(1-a)).*(t2.*Yp).^(1-a) - b.*(1./(2*Y)).*((t1.*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((t2.*Yp)-Fe2).^2)./...  
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fe2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fe1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fe2).^2)));

L_P=@(b) (1/(b*sigma*sqrt(2*pi)))*exp(-((log(b)-mu)^2)/(2*sigma^2)).*...
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*(t1.*Y).^(1-a)+Pfs*(1/(1-a)).*(t2.*Yp).^(1-a) - b.*(1./(2*Y)).*((t1.*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((t2.*Yp)-Fp2).^2)./...  
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fp2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fp1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fp2).^2)));

L_L=@(b) (1/(b*sigma*sqrt(2*pi)))*exp(-((log(b)-mu)^2)/(2*sigma^2)).*...
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*(t1.*Y).^(1-a)+Pfs*(1/(1-a)).*(t2.*Yp).^(1-a) - b.*(1./(2*Y)).*((t1.*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((t2.*Yp)-Fl2).^2)./...  
(exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((0/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((0/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((1/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((1/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((2/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((2/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((3/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((3/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((4/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((4/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((5/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((5/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((0/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((0/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((1/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((1/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((2/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((2/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((3/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((3/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((4/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((4/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((5/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((5/6)*Yp-Fl2).^2)...
 +exp(Pss*(1/(1-a)).*(Y).^(1-a)+Psf.*(1/(1-a)).*((6/6).*Y).^(1-a)+Pfs*(1/(1-a)).*((6/6).*Yp).^(1-a)- b.*(1./(2*Y)).*(((6/6)*Y)-Fl1).^2 - b.*(1./(2*Yp)).*((6/6)*Yp-Fl2).^2)));

L = [integral(L_E,0,inf,'ArrayValued',true), integral(L_P,0,inf,'ArrayValued',true), integral(L_L,0,inf,'ArrayValued',true)];
Li  = lambda_e*L(:,1) + lambda_l*L(:,3) + (1-lambda_l-lambda_e)*L(:,2);

likelihood = sum(log(Li))

FV_e_prob =  (lambda_e.*L(:,1))./Li;
FV_l_prob =  (lambda_l.*L(:,3))./Li;
FV_p_prob =  ((1-lambda_e-lambda_l).*L(:,2))./Li;
FV_prob=[L(:,1) L(:,2) L(:,3) FV_e_prob FV_p_prob FV_l_prob D(:,4)]; 

T_IN=table(FV_prob);
filename= 'IN_probdata_CRRA.xlsx';
writetable(T_IN,filename,'Sheet',1,'Range','D1')